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ABSTRACT 

Context. We study the physics of the solar tachocline (i.e. the thin transition layer between differential rotation in the convection zone and quasi 

■ uniform rotation in the radiative interior), and related MHD instabilities. 

' Aims. We have performed 3-D MHD simulations of the solar radiative interior to check whether a fossil magnetic field is able to prevent the 
[ spread of the tachocline. 

Methods. Starting with a purely poloidal magnetic field and a latitudinal shear meant to be imposed by the convection zone at the top of the 
radiation zone, we have investigated the interactions between magnetic fields, rotation and shear, using the spectral code ASH on massive 
parallel supercomputers. 

»C< Results. In all cases we have explored, the fossil field diffuses outward and ends up connecting with the convection zone, whose differential 
Oh' rotation is then imprinted at latitudes above « 40° throughout the radiative interior, according to Ferraro's law of isorotation. Rotation remains 
^ uniform in the lower latitude region which is contained within closed field lines. We find that the meridional flow cannot stop the inward 
' progression of the differential rotation. Further, we observe the development of non-axisymmetric magnetohydrodynamic instabilities, first 
| due to the initial poloidal configuration of the fossil field, and later to the toroidal field produced by shearing the poloidal field through the 
C$ , differential rotation. We do not find dynamo action as such in the radiative interior, since the mean poloidal field is not regenerated. But the 
^ instability persists during the whole evolution, while slowly decaying with the mean poloidal field; it is probably sustained by small departures 
. — from isorotation. 

■ Conclusions. According to our numerical simulations, a fossil magnetic field cannot prevent the radiative spread of the tachocline, and thus it 
S-h ' is unable to enforce uniform rotation in the radiation zone. Neither can the observed thinness of that layer be invoked as a proof for such an 

internal fossil magnetic field. 
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1. Introduction 

One of the great achievements of helioseismology was to map 
the internal rotation of the Sun (Brown et al. 1989). It revealed 
that the convection zone rotates differentially, with the angular 
velocity varying little in depth, and that the radiative interior 
is in quasi-uniform rotation. The transition between these two 
regimes occurs in a thin layer, the tachocline, which seems to 
straddle the boundary between the convection zone and the ra- 
diation zone. It is still a matter of debate whether the whole 
tachocline is turbulent, due to thermal convection, or only part 
of it. 

Here we shall assume that at least some differential rotation 
is applied on the top of the radiation zone proper, excluding the 
subadiabatic overshoot region. Then, in the absence of other 
physical processes, this differential rotation should spread into 
the radiation zone, as was shown by Spiegel & Zahn (1992). 

Send offprint requests to: A. S. Brun 



That spread is due to the thermal diffusion of the tempera- 
ture fluctuation which is generated in latitude, along isobars, 
by the differential rotation. These are tightly linked by the 
baroclinic relation; therefore, as the temperature fluctuation ex- 
pands downwards, so does the differential rotation. The process 
is not just a diffusion in this stably stratified medium, since it 
also involves meridional currents: the temporal variation of all 
fluctuations is proportional to the square of the Laplacian, and 
some refer to it as a hyper-diffusion. Viscosity too contributes 
to the spread, but only slightly. The whole process is very slow, 
but in the present Sun the tachocline should reach down to 0.3 
in radius, in the absence of any other physical process. 

Since this is not observed (Charbonneau et al. 1998; 
Corbard et al. 1999), one has to invoke a mechanism that tends 
to suppress the differential rotation in latitude, which is the 
cause of the spread. In Spiegel & Zahn's model, this is achieved 
by a turbulent anisotropic viscosity, with much stronger trans- 
port in the horizontal than in the vertical direction, which halts 
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the spread of the tachocline. Established originally in the thin 
layer limit, this property was confirmed by 2-dimensional sim- 
ulations performed by Elliott (1997). The fact that such tur- 
bulence acts to reduce the angular velocity gradient has been 
questioned by Gough & Mclntyre (1998), who cite geophysical 
examples where horizontal turbulence drives the system away 
from uniform rotation. On the other hand, in Couette-Taylor 
experiments turbulence reduces the angular velocity gradient, 
when angular momentum increases outward and the flow is lin- 
early stable; in other words it acts to suppress the shear which 
is responsible for the instability (Richard & Zahn 1999). This 
trend appears also in numerical 3-dimensional simulations of 
the turbulent overshoot region which have been carried out by 
Miesch (2003). 

Another way to oppose the spread of the tachocline is 
by magnetic stresses. The first to address the problem were 
Riidiger & Kitchatinov (1997). In their model, a steady mag- 
netic field confined in the radiation zone opposes the vis- 
cous spread of the differential rotation applied on the top; the 
tachocline reduces then to a very thin Ekman-Hartmann layer. 
MacGregor & Charbonneau (1999) confirmed that result, but 
they also considered another case, where the poloidal field 
was allowed to thread into the differentially rotating convection 
zone. A rather weak field then suffices to imprint the differential 
rotation throughout the radiation zone, thus enforcing Ferraro's 
law Bp • V£2 = 0, where B p is the mean poloidal magnetic field 
and Q. the angular velocity (Ferraro 1937). These models, how- 
ever, were rather crude, since they did not allow for the Ohmic 
diffusion of the poloidal field; moreover, they ignored thermal 
diffusion and the associated meridional circulation. 

The model outlined by Gough & Mclntyre (1998) includes 
most physical processes that are likely to operate in the upper 
radiation zone: the tachocline circulation is driven by thermal 
diffusion, and its spread is halted by the fossil field that resides 
beneath. The field does not penetrate into the tachocline, except 
perhaps in a narrow region of up-welling flow, and therefore 
it does not transmit the differential rotation of the convection 
zone to the deep interior. Actual calculations of this model have 
been carried out by Garaud (2000, 2002), directly seeking the 
stationary solution sketched by Gough & Mclntyre, although 
in her case the tachocline circulation is driven through Ekman 
pumping; some of her solutions show a tendency for mag- 
netic confinement, but a substantial differential rotation persists 
throughout the radiation zone, particularly at high latitude. 

The models discussed above favor a slow version of the 
tachocline with a typical ventilation time of the order of 10 6 yr. 
Gilman and collaborators (Gilman & Fox 1997; Dikpati, Cally 
& Gilman 2004 and references therein) developed a series of 
models that showed that the tachocline could become unstable 
through magnetic instability of toroidal structures embedded 
in it, resulting in a fast latitudinal angular momentum trans- 
port that suppresses the shear and limits its inward diffusion. 
Forgacs-Dajka and Petrovay (2002) considered the effect of an 
oscillating dynamo field of 22 yr period and found that it too 
could suppress the spread of the solar tachocline. Both models 
favor a rather fast tachocline with ventilation times of the order 
of a rotation period or the duration of the solar cycle (Gilman 
2000). 



In the present paper we do not deal with such short 
timescales, and we examine a model of the solar radiative zone 
and tachocline similar to that of Gough & Mclntyre. But we 
consider the time-dependent problem, starting from various ini- 
tial conditions, for two reasons: i) both the tachocline spread 
and the diffusion of the magnetic field proceed very slowly, 
and thus it is not clear that a stationary solution can be reached 
by the age of the Sun; ii) such a stationary solution may not 
be unique, but may actually depend on the initial conditions. A 
similar approach has been taken by Sule et al. (2005); we differ 
from their work in that we let the poloidal field diffuse, and in 
our case the tachocline circulation is driven mainly by thermal 
diffusion, as in the Sun. 

Moreover, we use for our simulations the 3-dimensional 
code ASH (see further down) and resolve the Alfven cross- 
ing time; this enables us to describe non-axisymmetric MHD 
instabilities which could lead to drastic reconfigurations of the 
magnetic field, as was recently demonstrated numerically by 
Braithwaite and Spruit (2004; see also Braithwaite & Nordlund 
2005), following the pioneering work of Tayler (Tayler 1973; 
Pitts & Tayler 1985). 

2. The model 

We make use of the well-tested hydrodynamic ASH code 
(Anelastic Spherical Harmonic; see Clune et al. 1999; Miesch 
et al. 2000; Brun & Toomre 2002) which was originally de- 
signed to model the solar convection zone. It has since been 
extended to include the magnetic induction equation and the 
feedback of the field on the flow via Lorentz forces and Ohmic 
heating (see Brun, Miesch & Toomre 2004 for more details). 
Here we apply it to examine the nonlinear interaction between 
a shearing flow, i.e. the differential rotation imposed by the con- 
vection zone on the top of the radiative interior, and a fossil 
magnetic field. Our goal is to assess whether such a field can 
hinder the spread of the tachocline deep into the radiation zone 
and thus keep it thin, as was revealed by helioseismic inver- 
sions (Charbonneau et al. 1998; Corbard et al. 1999). 

Our numerical model is a simplified portrayal of the so- 
lar radiation zone: solar values are taken for the heat flux, 
rotation rate, mass and radius, and a perfect gas is assumed. 
The computational domain extends from r DO t = 0.35 R e to 
lop = R = 0.72 R Q (Rq is the solar radius); we thus focus on the 
bulk of the stably stratified zone, excluding the nuclear central 
region, and we ignore its possible back reaction on the convec- 
tive envelope. The vertical extent of the computational domain 
is therefore D = 2.5 x 10 10 cm and the background density 
varies across the shell by about a factor of 17. 

2.1. Governing equations 

The ASH code solves the full set of 3-D MHD anelastic equa- 
tions of motion (Glatzmaier 1984) in a rotating spherical shell 
on massively-parallel computer architectures (Brun, Miesch 
& Toomre 2004). Depending on the sign of the initial mean 
entropy gradient considered in the model, either negative or 
positive, convection or radiation zones can be simulated. The 
anelastic approximation is used to retain the important effects 
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of density stratification without having to track the sound waves 
(i.e. dp/dt = 0). The MHD anelastic equations are fully non- 
linear in velocity and magnetic field variables, but under the 
anelastic approximation the thermodynamic variables are lin- 
earized with respect to a spherically symmetric and evolving 
mean state, with density p(r,t), pressure P(r,t), temperature 
T (r, t) and specific entropy S (r, t). Fluctuations about this mean 
state are denoted by p, P, T, and S . The resulting equations are: 



impose a differential rotation meant to be enforced by the 
turbulent convection zone. Following helioseismic inversions 
(Thompson et al. 2003), we apply the following law: 



Q(r top , 0) = A + B cos 2 6 + C cos 4 G. 



(8) 



The uniform rotation which has the same specific angular mo- 
mentum (Gilman et al. 1989) is 



V-(pv) = 0, 
V B = 0, 



(1) 
(2) 



1 3 
Oo=A + -B+ —C; 
5 35 



(9) 



pi — +(v-V)v + 2Q 



x vj = - 



VP + pg+ — (V x B) x B 

4n 



V £)- [VP-pg], 



(3) 



pT— +pTv-V(S +S) = V • [pcJRVT + kVT)] 
at 



we shall apply this uniform rotation initially to the whole radia- 
tion zone, and therefore there will be no net viscous flux across 
the boundary. 

For the remaining top and bottom boundary conditions, we 
impose: 



+ +2pv [ej/e,j - 1 /3(V ■ v) 2 J 1 ■ Impenetrable walls at the top and bottom: 



4717/ >2 __ 



(4) 



~dt 



= Vx(vxB)- V x (t/VxB), (5) 



where v = (v r , vg, Vq) is the local velocity in spherical coordi- 
nates in the frame rotating at constant angular velocity £2 , g 
is the gravitational acceleration, B = (B r , B g , B^) the magnetic 
field, j = c/4tt(V x B) the current density, and c p the specific 
heat at constant pressure. D is the viscous stress tensor, involv- 
ing the components 

2> y = -20v[*y-l/3(V-v)*y], (6) 

where is the strain rate tensor. 

The radiative diffusivity k applied to the mean temperature 
gradient Vf takes the solar value, hence the radiation flux is 
that of the Sun to keep the model consistent. But the diffusion 
coefficient k operating on the temperature fluctuations will be 
enhanced in our simulations, as explained below, and so also 
will the viscosity v and the Ohmic diffusivity 77. A volume heat- 
ing term pe is also included in these equations for completeness 
(see Brun, Browning & Toomre 2005). 

Finally we use the linearized equation of state, assuming 
the ideal gas law: 

p P T yP c p 

where y is the adiabatic exponent. We neglect the composi- 
tion gradient due to element settling. The reference or mean 
state, which is indicated by overbars, is derived from a one- 
dimensional solar structure model (Brun et al. 2002); it is fre- 
quently updated with the spherically-symmetric components 
of the thermodynamic fluctuations as the simulation proceeds. 
The model begins in hydrostatic balance so the bracketed 
term on the right-hand-side of equation (3) initially vanishes. 
Departures from this hydrostatic balance are found to be small. 

2.2. Integration domain, boundary conditions 

The system of equations to be solved with ASH requires 12 
boundary conditions. At the top of the numerical domain we 



V r — 0| r=rboti r top , 

2. Stress free at the bottom: 

!(?)-!(?)-*- 

3. Constant entropy gradient at the top and bottom: 
dS_ _ 

4. We match the magnetic field to an external potential field 
at the top and bottom: 

B = V(D -» AO = 0U bo „ r , op . 

Our choice of magnetic boundary conditions ensures that no 
magnetic torques are applied to the numerical domain. 

2.3. Choosing the physical parameters 

In Gough & Mclntyre's model, a magnetic boundary layer, 
or magnetopause, separates the magnetic interior from the 
tachocline void of magnetic field; its thickness 6 is determined 
by the relative strength of the two diffusion processes and that 
of the magnetic field: 



\R) L 4 



2 (f A ) 2 



(10) 



where R = r top , and L a 4.5 is related to the latitudinal structure 
of the tachocline. Here we have introduced the Alfven crossing 
time fA = R/Va = R y/4np/B; f B = R 2 h is the magnetic diffu- 
sion time, and ?es = (N/£lo) (R Ik) the local Eddington-Sweet 
time. B is the magnetic field strength and N the buoyancy fre- 
quency. 

Expressing that the diffusion of the magnetic field is halted 
by the downflow of the tachocline circulation, whose ventila- 
tion time is given by 



(11) 
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Fig. 1. Radial dependence of the Prandtl P r , magnetic Prandtl 
P m and Roberts Q numbers in the Sun based on the microscopic 
values of the diffusivities for a ionized hydrogen gas (Zeldovich 
et al. 1983). 



Parameter 


Symbol 


Sun 


Simulation 


mean thermal diffusivity 


K 


i n7 
10 


1 n7 
10 


thermal diffusivity 


K 


1 n7 
10' 


8 10 12 


magnetic diffusivity 


n 




8 10 lu 


viscosity 


V 


30 


8 10 9 


Prandtl number 


v/k 


3 10" 6 


10 3 


Roberts number 


j]/k 


10 


i n-2 
10 


magnetic Prandtl nb 


v/i] 


310~ 2 


io- 1 


Froude number 


(Ho/AO 2 


21Q- 6 


io~ 4 


Ekman number 


v/2Q R 2 


210~ 15 


610- 7 


rotation frequency 


Sh 


3 10 6 


310- 6 


Alfven time (for 1 G) 


t A 


5 10 10 


5 10 10 


Ohmic diffusion time 


t B 


2.5 10 18 


3.125 10 10 


viscous diffusion time 


ty 


8.333 10 19 


3.125 10" 


Eddington- Sweet time 


tss 


1.25 1 20 


3.125 10 12 



Table 1. Typical values of the relevant parameters in the upper 
radiation zone of the Sun, and values adopted for the model 
at r top = O.72/? (in cgs units, when applies). In the simula- 
tion, the thermal diffusivity k operating on the mean tempera- 
ture gradient takes its solar value; only the diffusivity k applied 
on the temperature fluctuations is increased. 



one obtains the following relation between the thickness 6 of 
the magnetopause and A, that of the tachocline : 

\R) L\Rl ^{Qo)' 

where Q/Qo ~ 1/5 is a measure of the differential rotation. 

We note that the size of these boundary layers does not de- 
pend too sensitively on the value of the diffusivities and on 
the strength of the magnetic field, which means that we will 
not alter too much the problem when we increase these values 
to meet the numerical requirements. Figure ^shows the vari- 
ation with depth of the non-dimensional numbers that charac- 
terize the microscopic diffusivities in the Sun, and typical val- 
ues for the relevant parameters in the upper radiation zone are 
given in Table [0 Due to limitations in computing resources, 
no simulation achievable now or in the near future can hope 
to directly use microscopic values for the diffusivities. In our 
model we had to increase all diffusivities to be able to resolve 
the smallest scales (i.e. the thickness of the Ekman layer, and 
the size of the turbulent eddies induced by the instabilities de- 
scribed later in ©, but we took care to respect their hierar- 
chy. We chose the viscosity v as small as allowed by our nu- 
merical resolution, in order to ensure that the tachocline be 
driven by thermal diffusion, as in the Sun, rather than through 
Ekman pumping. Another way to look at the problem is to 
make sure that the spread due to thermal hyper-diffusion, which 
evolves as A/R w (/Aes) 1 ^ 4 , still exceeds the viscous spread 
A/R a (t/t v ) 1/2 at the age of the Sun (with t v = R 2 /v), which 
measured in Eddington-Sweet time is £ /fEs = 1.2 10~ 3 ; this 
is guaranteed if the Prandtl number is less than v/k = 3 10~ 3 , 
hence our choice of v/k = 10~ 3 . 

The Froude number Q.q/N was also increased to meet the 
numerical requirements; since this enhances the ventilation 
speed, it may also account for the fact that the Sun was ro- 
tating more rapidly in the past. With these parameters, taking 



a field of lkG yields a tractable thickness for the magnetic 
boundary layer: 5/R = 2.2 10~ 2 according to (II 0i . From fl!2i 
the corresponding tachocline thickness is then A/R = 2.1 10~ 2 , 
and the magnetic Reynolds number based on the Alfven speed 
R m = ts/tA will reach a few 10 3 in the computational domain. 

2.4. Numerical resolution 

The global code ASH solves in a 3-D rotating spherical shell 
the set of equations (1-5) in the anelastic approximation. This 
has the great advantage that the Courant stability criterion is 
based on the slow motions present in the radiative zone rather 
than on the sound speed, thus allowing for much larger time 
steps. In the MHD context, the anelastic approximation filters 
out fast magneto-acoustic waves but retains the Alfven and 
slow magneto-acoustic modes. 

The ASH code is based on the so-called pseudo-spectral 
numerical method to solve the system of partial differen- 
tial equations introduced above (Glatzmaier 1984; Canuto et 
al. 1995; Boyd 1989; Clune et al. 1999). In this numerical 
method, the main variables are projected onto orthogonal func- 
tions/polynomials rather than discretized on a mesh grid as 
with finite differences. The advantage of this method is that in 
spectral space, derivatives are straightforward multiplications, 
rather than variations as in a finite difference scheme, which 
are numerically less precise. Since we want to solve the MHD 
equations in spherical geometry, it is natural to use the spheri- 
cal system of coordinates (r, 6, <p) and to project the main vari- 
ables on spherical harmonics Y( m {9, (/>). This approach has the 
advantage that the spatial resolution is uniform everywhere on 
a sphere when a complete set of spherical harmonics is used up 
to some maximum in degree t (retaining all azimuthal orders 
m < £ in what is known as triangular truncation). Gaussian 
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quadrature imposes that we use the relation i max — (2Ng - l)/3, 
and N<p = 2Ng for the horizontal resolution (Boyd 1989). For 
the radial direction we choose to project our functions onto 
Chebyshev polynomials T n (r), using N r Gauss-Lobatto colo- 
cation points (Canuto et al. 1995). Typical grid resolution for 
our models is N r = 193, Ng = 128 and Nq = 256; the time-step 
was chosen such as to resolve the gravity modes. 

In order to ensure that the mass flux and the magnetic field 
remain divergence-less to machine precision throughout the 
simulation, we use a toroidal-poloidal decomposition: 



pv = V x V x (We r ) + V x (Ze,.), 
B = V x V x (Ce r ) + V x (Ae r ) . 



(13) 
(14) 



More details on the numerical method can be found in Clune et 
al. (1999) and Brun, Miesch & Toomre (2004). 

2.5. Initial conditions 

At t = the computational domain is motionless, except for a 
uniform rotation Q = £2o defined in l|9}, and for the differential 
rotation imposed at the top boundary r = r top , where the fol- 
lowing values have been chosen for A, B and C in expression 
©: A/2n = 456, B/2n = -42 and C/2n = -72 nHz. 

Initially we apply a purely axisymmetric meridional field 



B r e, + BgCg, with 



go ff¥ 
r 2 sin 8 38 ' 



Bg = - 



go gg. 

r sin 8 3r ' 



(15) 



*¥(r, 8) is constant on field lines, and is related to the azimuthal 
component of the vector potential: *P = r sin dA^. For a dipolar 
field *P oc sin 2 8R 3 jr; here we take instead 



¥ = (r/R) z (R B - rf sin z 8 for r < R B 
¥ = for r > R B , 



(16) 



which confines the field in the sphere of radius Rb- 1 

This parameter Rb thus determines the degree of confine- 
ment of the initial field; we shall consider the 3 cases: 

- I. R B » r top - the field threads into the convection zone, 

- II. Rb = r top = O.72R - the field is just barely confined in 
the radiation zone, 

- III. Rb - O.64/? - the initial field is buried deep in the 
radiation zone. 



3. Results 

The interaction of rotation and magnetic field in the solar ra- 
diative zone is found to be very sensitive to the initial magnetic 
field configuration, especially in the early phases of the tempo- 
ral evolution; we shall examine in turn the cases I, II and III de- 
fined above. To ease the comparison with the real Sun, we have 
rescaled the time by the ratio fEs(Sun)/fEs(simulation) = 4 10 7 ; 
thus from here on we quote the time in that 'solar equivalent 
unit'. 



1 Charbonneau and MacGregor (1993) took a similar form for *P; 
we differ from them by enforcing the continuity of B B at R B , thus sup- 
pressing there the current sheet. 



3.1. Case I: Initial magnetic field threading into the 
convection zone 

In this first case, at t — the magnetic field lines are already 
threading into the differentially rotating convection zone, as 
shown in Fig. |2^. The torque applied by the shear at the top 
of the domain is quickly transmitted throughout the radiation 
zone, first at high latitudes where the direction of the field lines 
is almost radial (Fig. |2j)) and then more slowly at lower lati- 
tudes (Fig. 2J;). Such a strong differential rotation in latitude 
in the solar radiative interior is not observed in helioseismic in- 
versions (Thompson et al. 2003; Couvidat et al. 2003). Clearly, 
in this case the magnetic stresses do not oppose the inward 
propagation of the shear but they favor it. (Similar results have 
been found in 2-D by MacGregor & Charbonneau 1999.) The 
transport of angular momentum by Maxwell stresses is very 
fast, and prevents the magnetopause from forming. Meridional 
flows do exist but they are not found to play a significant role: 
they are unable to oppose the spread of the latitudinal shear. 
Their amplitude in the simulation is about 1 m/s, which cor- 
responds (after rescaling) to a real solar ventilation speed of 
10~ 7 m/s, in excellent agreement with the theoretical predic- 
tion (Gough & Mclntyre 1998). 

The interaction between the poloidal field and the shear 
leads initially to a fast increase of the axisymmetric toroidal 
field, which then saturates as the isocontours of Q. become more 
and more aligned with the axisymmetric poloidal field, tending 
to Ferraro's law of isorotation. Its strength, over a short interval 
of time, is larger than that of the axisymmetric poloidal field, 
before equilibrating at about the same amplitude. These mean 
toroidal fields are found to become unstable to a m — 1 insta- 
bility. Deeper down and in the early phase of evolution, we ob- 
serve the fast grow of high-w non-axisymmetric perturbations, 
which is due to a MHD instability of the purely poloidal field 
configuration considered. Both the low and high-ra instabilities 
will be described in more detail in |4] 

3.2. Case //: Initial field barely confined in the radiation 
zone 

In the second case, the magnetic field lines do not penetrate 
into the convection zone but they come very close to it, as can 
be seen in Fig. [3^. At the beginning of the simulation, the vis- 
cous torque associated with the imposed shear at the top of the 
domain does not drag the magnetic field lines, and as a conse- 
quence no physical processes are available to quickly transmit 
the latitudinal shear. This particular configuration of the mag- 
netic field has been chosen in the 2-D calculations of Riidiger 
& Kitchatinov (1997) and MacGregor & Charbonneau (1999) 
to demonstrate the confinement of the tachocline. If the initial 
poloidal magnetic field is kept fixed in time, as in their calcu- 
lations, its peculiar topology is quite efficient at stopping the 
inward diffusion of the latitudinal shear. However, in our simu- 
lations the magnetic field is allowed to diffuse outward, where 
it connects with the imposed rotation (Fig.[3j)). As the more in- 
tense field lines reach the top boundary, and slowly drift toward 
the polar regions while becoming more radial, isorotation is es- 
tablished here also, first at high latitude and then at lower lati- 
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Fig. 2. Case I. Temporal evolution, over 860 Myr in equivalent solar units, of the angular velocity Q. (color contours and white 
lines) and the mean axisymmetric poloidal field (superimposed black lines, whose labels refer to the field strength in the equatorial 
plane). The time is measured in solar equivalent units, i.e. rescaled by the ratio fes(Sun)/fEs(simulation) = 4 10 7 . Initially, the 
field already penetrates into the convection zone; it quickly imprints the differential rotation imposed at the top on the whole 
radiation zone. 




Fig. 3. Case II. Temporal evolution, over 950 Myr in solar equivalent units, of the angular velocity CI (color contours and white 
lines) and the mean axisymmetric poloidal field (black lines). The initial field does not penetrate into the convection zone, but 
the outcome is the same as in case I. 



tude (Fig.|5J;). It takes about 8 times longer for Ferraro's law to 
be achieved in case II than in case I. However this is still a fast 
process, accomplished at an age of 0.2 Gyr in equivalent solar 
units. Here again the internal magnetic field fails to prevent the 
inward spread of the shear because the magnetopause does not 
have the time to form, and as a consequence there is no merid- 
ional circulation to oppose the spread of the latitudinal shear by 



horizontally bending the magnetic field lines. Meridional flows 
do exist but here again they do not play a significant role. 

As soon as the field lines make contact with the shear, the 
situation we observed in the first case is again realized, but the 
toroidal field takes about 8 times longer to reach the same level 
as in case I. The energy in the mean toroidal field saturates as 
the isocontours of Q become more and more aligned with the 
axisymmetric poloidal field. As in case I, we find that the mean 
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Fig. 4. Case III. Temporal evolution of the angular velocity £2 (color contours and white lines) and the mean axisymmetric 
poloidal field (superimposed black lines), which initially is buried below the convection zone, a-d) sequence spanning 4.7 Gyr 
(in solar equivalent units) of the dipolar magnetic field in the presence of rotation and shear. The field lines gradually connect with 
the imposed top shear, thereby enforcing differential rotation at all depths above latitudes greater than « 40°, at the (equivalent) 
age of the Sun. Note also in frames c-d the distortion near the poles of the lines of constant O, which is associated with the Tayler 
instability. 




Fig. 5. a-d) Case III. Temporal evolution of the associated toroidal field (color contours) and of the meridional circulation 
(superimposed black lines) at t= 81, 280 Myr, 1.3 and 4.7 Gyr (solar equivalent units). We clearly see the early formation of a 
magnetopause with an intense mean (axisymmetric) toroidal field, and then the development at high latitudes of a toroidal field 
associated with the progression of the latitudinal shear into the radiation zone. Note the highly distorted contours in c-d, a sign 
that the Tayler instability has reached high amplitude, and the strongly radially squeezed counter-cells near the surface at low 
latitudes. 

toroidal field becomes unstable to a m = 1 instability at high that case II leads to the same state of differentially rotating ra- 
latitude, and that a rapid non-axisymmetric MHD instability of diative interior as case I, except that this occurs later. Hence 
the mean poloidal field develops in the early phase. The tim- we do not confirm the results of Riidiger & Kitchatinov (1997) 
ing of these instabilities varies between the 3 cases, being the and MacGregor & Charbonneau (1999), that a barely confined 
strongest for case III (see below §4). We can thus conclude field is able to prevent the spread of the tachocline. The tempo- 
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ral evolution of the mean poloidal field plays a crucial role here 
by allowing the field to connect to the imposed upper shear. 



3.3. Case III. Initial field buried deeply in the radiation 
zone 

In this last case, the magnetic field lines are initially confined 
below r - Rb - O.64R , to let the tachocline penetrate some- 
what into the interior before encountering the fossil field (Fig. 
|4^). When the field lines make contact with the shear (see 
Fig-EbX we notice a fast increase of the mean toroidal energy 
(TME) in a thin band extending over a large range in latitude, 
which corresponds to the magnetopause anticipated by Gough 
& Mclntyre (see Fig. [5^). In case III, it takes longer for the 
TME to reach a level comparable to that of cases I and II, as 
expected since at the beginning the mean poloidal field does 
not interact with the shear. 

At 280 Myr, the meridional circulation in the tachocline 
consists of three superposed counter-rotating cells in each 
hemisphere, of which the largest spans all latitudes, like a 
dipole. The pattern resembles that obtained by Elliott (1997) 
in the non-magnetic case. Here also the circulation is octopo- 
lar close to the top boundary, with a tilted up-welling located 
around the latitude of 35°. This meridional circulation is un- 
able to prevent the outward diffusion of the magnetic field lines, 
even in the down-dwelling regions (Figs. 03-d); it does not suc- 
ceed in squeezing the field lines in narrow up-wells so as to im- 
print below only a small amount of differential rotation. The di- 
rect consequence of this inefficiency is that the field lines even- 
tually connect to a large portion of the shear imposed at the top 
of the domain by diffusing through the magnetopause. This re- 
sults in an upward shift of the magnetopause, which is the site 
of a strong generation of toroidal field (via the Q-effect) - notice 
the migration of B<f, between Figs. [5^ and b. Then as the field 
lines gradually connect to the shear imposed at higher latitudes, 
becoming more radially oriented, they communicate the differ- 
ential rotation of the convection zone to all depths via Maxwell 
stresses, as in the two previous cases. At 4.7 Gyr, the merid- 
ional circulation has shrunk into two very thin dipolar cells 
in each hemisphere, one being hardly distinguishable near the 
equator, and the differential rotation has invaded all latitudes 
above 40° from top to bottom. 

We have run Case III with different values of the confine- 
ment depth Rb and of the Prandtl number v/k, which allows us 
to conclude that the speed at which the differential rotation pen- 
etrates into the radiative interior depends mainly on the time it 
takes for the poloidal field lines to connect with the convec- 
tion zone: the deeper the field is buried initially and the longer 
it can maintain uniform rotation within its closed field lines, 
while allowing for a thicker tachocline. Thus Ferraro's law of 
isorotation, with a radiative interior rotating differentially, is 
the most likely outcome of the interaction of the poloidal field 
with the latitudinal shear. The initial state of solid rotation only 
lasts as long as the fields lines have not diffused outward and 
connected with the convection zone. 




2X10 



3x10 



4X10 



Fig. 6. Case III: temporal evolution of kinetic and magnetic 
energies in the presence of rotation and shear. We clearly see 
the 2 phases of non-axisymmetric instability (characterized by 
the fast increase of FKE and FME): the first, for t < 1 Gyr, is 
associated with the unstable configuration of the purely dipolar 
initial field, and the second, for t > 1 Gyr, with the toroidal 
field produced by winding up the poloidal field through the dif- 
ferential rotation. The time is given in solar equivalent units, i.e. 
rescaled by the ratio fEs(Sun)/fEs(simulation). Diamonds have 
been superimposed on the FME curve to indicate the temporal 
sampling used in Figs.0and|8] 



4. Instabilities 

Our three-dimensional simulations reveal another interesting 
and fundamental aspect of the problem, namely the occurrence 
of non-axisymmetric MHD instabilities in the radiation zone. 
These have been studied in detail through linear analysis by 
Tayler and collaborators (Markey & Tayler 1973, 1974; Tayler 
1980; Pitts & Tayler 1985) and by Wright (1973); an excellent 
review of the subject has been written by Spruit (1999). Their 
main finding was that purely poloidal or purely toroidal fields 
are unstable to non-axisymmetric perturbations, and that sta- 
ble configurations must therefore be of mixed poloidal/toroidal 
type, with comparable field strength. Since then, such config- 
urations have indeed been identified in recent numerical simu- 
lations performed by Braithwaite & Spruit (2004, 2006) and 
Braithwaite & Nordlund (2006). Although our primary goal 
here was to examine the possibility that the solar tachocline 
could be confined by a fossil magnetic field, we feel that our 
calculations have also a broader scope, and that they may con- 
tribute to better understanding the behavior of stellar magnetic 
fields. We shall concentrate here on case III; cases I and II show 
similar behavior. 

A useful quantitative diagnostic is provided by Fig. [6] 
which describes the temporal evolution of various components 
of the kinetic (KE) and magnetic (ME) energies, averaged over 
the whole computational domain. PKE & PME designate re- 
spectively the mean (axisymmetric) poloidal components of 
KE and ME, TKE & TME their mean toroidal components, and 
FKE & FME their non-axisymmetric components (see Brun et 
al. 2004 for their analytic expressions). In Fig.0we display the 
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Fig. 7. Case III. Temporal evolution of the azimuthal magnetic field (from which we have substracted the m = component), 
shown in horizontal Mollweide projection as a function of depth (respectively from left to right: r = 0.7, 0.55 and 0.37^ o ), at 
various times in the evolution as indicated in solar equivalent units. First, the high-m instability of the poloidal field appears near 
the equator at the bottom of the domain, and later the low-m instability of the toroidal field occurs at all depths near the poles. 
Note the remarkably clean m = 1 structure at the top in the frames for 0.85 and 1.3 Gyr. Between ~ 500 Myr and 850 Myr there 
is a clear coexsistence of both instabilities in the simulation. 
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maps, in Mollweide projection, of the non-axisymmetric part 
of the toroidal field, at different depths and epochs; they will 
serve to locate and to characterize the instabilities. 

^From the onset, starting from an extremely weak non- 
axisymmetric velocity field, FKE & FME rise exponentially 
by many orders of magnitude at a fast pace (the e-folding time 
is w 10 Myr), and then they saturate around 700 Myr at a level 
1000 times below the energy of the mean poloidal field (PME). 
As can be seen in Fig. [7] the fluctuations are located at low lat- 
itude and at the bottom of the domain; we conclude that it is 
the instability described by the pioneers quoted above, which 
occurs when the field is purely poloidal. Figure|S]displays the 
energy spectrum in m of the fluctuations of the azimuthal field 
B$ at the bottom of the domain: one sees clearly that the max- 
imum growth occurs at rather high azimuthal wavenumber: 
m as 40. This is in qualitative agreement with Wright (1973) 
and Markey & Tayler (1974), who found that the growth-rate 
should steadily increase with m in the non-dissipative case. 
Note in Fig. [5] that no mean toroidal field is produced in the 
unstable region, presumably because the magnetic helicity was 
initially zero. 



Bphi Spectra at r=263.30 Mm 
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Fig. 8. Distribution over azimuthal wave-number m (and 
summed of all I > m) of the magnetic energy in B$, close to 
the bottom of the domain. The two instability phases have dis- 
tinct signatures: the first, due to the purely poloidal field con- 
figuration, occurs at m > 20 with a maximum around m = 40, 
whereas the second, caused by the toroidal field, starts around 
450 Myr, and involves low m wave-numbers. At about 1.3 Gyr, 
the highest wave-numbers are eroded through Ohmic diffusion. 

In the meanwhile, two bands of toroidal field of opposite 
polarity appeared in the upper part of the domain, one at mid 
latitude and the other close to the rotation axis, with opposite 
sign in each hemisphere. This field is generated by the shear- 
ing of the poloidal field, which very soon encountered the dif- 
ferential rotation spreading down from the convection zone. At 
about 1 .2 Gyr the energy TME, averaged over the whole do- 
main, matches that of the mean poloidal field, meaning that lo- 
cally the toroidal field can be much stronger than the poloidal 



field, since it occupies a smaller volume. That toroidal field is 
unstable to low m perturbations, as can be seen both in Fig.Q 
and in the spectra of Fig. [3] Note that the mean toroidal field 
and the associated instability draw their energy from the differ- 
ential rotation, as indicated by the conspicuous inflection of the 
TKE curve in Fig. [5] at 1.1 Gyr. 

Thereafter all energies but one slowly decline, at a rate 
which is controlled mainly by the Ohmic dissipation of the 
mean poloidal field (PME): in our rescaled solar units, the e- 
folding time would be R 2 /2n 2 ri « 2 Gyr for a uniform mo- 
tionless sphere (Cowling 1957), which is compatible with what 
we see here. The decay of that field is particularly smooth, 
and there is no sign that it is regenerated or even affected by 
the instability. The mean toroidal field accompanies closely 
the decline of the mean poloidal field, and so do the non- 
axisymmetric components FKE and FME, with the kinetic en- 
ergy in these perturbations being about twice their magnetic 
energy. The only exception to that overall decline is the mean 
toroidal kinetic energy TKE, which steadily increases as the 
differential rotation spreads deeper and to lower latitudes. The 
kinetic energy of the meridional flow (PKE) is the smallest of 
all energies; it is mostly concentrated toward the top of the do- 
main and it remains almost constant during the whole evolu- 
tion. 

Recently Braithwaite & Spruit (2004) studied Tayler's in- 
stability in a polytropic star, by means of 3-D numerical simu- 
lations. They found that a random initial magnetic field relaxes 
in a few Alfven crossing times into a mixed poloidal/toroidal 
topology, by exciting a non-axisymmetric small scale instabil- 
ity. This stabilized field then slowly decays and migrates to- 
ward the surface. Our simulations differ here from Braithwaite 
& Spruit in that we include rotation, initially impose a large 
scale poloidal field, thus with zero helicity, and that we ap- 
ply differential rotation at the top the domain, which shears 
the slowly decaying poloidal field to produce a toroidal field, 
antisymmetrical with respect to the equator. Near the poles, 
our field is of the kind considered by Tayler (1973) and it is 
known to give rise to non-axisymmetric instabilities at low az- 
imuthal wave-number m; these are observed here in a fully de- 
veloped regime, over a broad spectrum of scales, suggesting 
that they could mix the chemical elements in the radiation zone. 
However, none of the instabilities observed in our simulations 
seems able to regenerate the mean (m = 0) poloidal field from 
which we started the simulations. This indicates that our sim- 
ulations cannot be considered as dynamos, since the dynamo 
loop (Bp — > B, — > Bp) is not fulfilled. We stress however that 
this model does not apply to the Sun, since such a fossil field 
alone does not allow for a thin tachocline, as we have seen. 

5. Discussion and conclusions 

We undertook these 3-D simulations primarily to check 
whether a fossil magnetic field can prevent the radiative spread 
of the solar tachocline, as was advocated by Gough and 
Mclntyre (1998). In all cases we have explored, the outcome 
is always that the poloidal field connects with the convection 
zone, and that it then imprints the differential rotation of that 
zone throughout the radiative interior, according to Ferraro's 
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law, at latitudes higher than about 40° (see Figs. 4c-d). When 
the initial poloidal field threads into the convection zone, the 
process takes only a few Alfven crossing times, as expected. 
When the field is confined in the radiation zone, the time re- 
quired for this depends on how deeply it is buried. Nearly uni- 
form rotation can only persist inside the outermost closed field 
line, hence at low latitude. Of course, a deeply buried field will 
not reach the convection zone by the age of the Sun, but then 
the tachocline will also have penetrated very far. 

The way the topology of the poloidal field enforces 
a given rotation profile was elucidated by MacGregor and 
Charbonneau (1999), but they did not let the field diffuse into 
the convection zone. The stationary solutions built by Garaud 
(2000) resemble ours, especially in what she calls the low field 
case which applies to the field strength chosen here; the main 
difference to our model is probably that her meridional circula- 
tion is driven by Ekman pumping. 

We believe that our simulations have captured the essence 
of the real problem, since they respect the similarity of time- 
scales characterizing the two competing processes, namely 
the Ohmic diffusion time of the mean poloidal field and the 
Eddington-Sweet time which rules the ventilation and the 
spread of the tachocline. In Gough & Mclntyre's model, an 
important role is played by the magnetopause, namely that the 
boundary layer that separates the poloidal field in the uniformly 
rotating interior from the tachocline void of magnetic field: it 
is the seat of the production of the toroidal field. According to 
JlOi its thickness 5 scales as B p 1 ^, and to allow its existence 
in spite of the enhanced diffusivities, i.e. in order to ensure 
8 <K r top , we had to increase the initial poloidal field B p by 3 or- 
ders of magnitude with respect to what would be required in the 
Sun. This could still be insufficient, since in the cases discussed 
here the two boundary layers have similar thickness, whereas 
in Gough & Mclntyre's model the magnetopause is much thin- 
ner than the tachocline. To check this point we ran a simulation 
with enhanced thermal diffusivity and reduced Ohmic diffusiv- 
ity, to increase the contrast A/5 by a factor of 3: the outcome 
was the same, and the differences barely noticeable, except for 
slightly less regular poloidal field lines. 

A more delicate point is whether the poloidal field lines ac- 
tually penetrate into the convection zone, as assumed in most 
models, including ours, or whether they are deflected by the 
turbulent motions, which are probably structured in plumes. 
Realistic simulations of penetrative convection are required to 
see how a magnetic field, anchored in the deep interior, will be- 
have there; we intend to tackle this problem in the near future. 

Since we used a three-dimensional code, we were able to 
observe the non-axisymmetric instabilities associated with our 
field configurations. In the early phase, they occur in a narrow 
equatorial belt at the base of the computational domain, and 
they are due to the purely poloidal field we have imposed ini- 
tially. Later on, they are present at all depths near the poles, 
where the shear of the differential rotation keeps generating a 
strong toroidal field. Their properties agree with the predictions 
of Tayler and his collaborators, which were recently reviewed 
and completed by Spruit (1999, 2002). 

An important result of our simulations is that these instabil- 
ities do not interfere with the mean poloidal field; they are able 



to distort the 'isogyres' (surfaces of constant angular velocity 
£2), thus affecting somewhat the production of the toroidal field, 
but even there they seem to have limited impact. Further we do 
not see a dynamo process occuring in our simulated radiative 
interior since the mean poloidal field is not regenerated and de- 
cays away. 

Based on their model, and on the observed thinness of the 
tachocline, Gough and Mclntyre (1998) concluded about the 
inevitability of a magnetic field in the Sun's radiative zone. We 
find here - but note the caveats above - that such a field seems 
unable to prevent the spread of that shear layer, and therefore 
that the thinness of the tachocline cannot prove or disprove the 
existence of a fossil field. This conclusion holds even when the 
aforementioned instabilities are taken into account, as in the 
new model presented by Mclntyre (2006), since these instabil- 
ities do not affect the mean poloidal field. 

Therefore other processes must be invoked to account for 
the thinness of the tachocline. Or one has to explain why there 
is no need for such a layer at the top of the radiation zone, and 
thus how uniform rotation is accomplished already at the base 
of the convection zone. A first step in that direction has been 
taken by Forgacs-Dajka and Petrovay (2000, 2002): they intro- 
duce the oscillating dynamo field in the region of convective 
penetration, and they succeed in smoothing out the differential 
rotation. 

But other possibilities are being actively explored. Talon 
and Charbonnel (2005; see also Charbonnel & Talon 2005) 
have shown that such uniform rotation can be achieved with in- 
ternal gravity waves emitted at the base of the convection zone, 
combined with meridional circulation and shear-induced turbu- 
lence in the bulk of the radiation zone. It is clear that in order to 
disentangle the theoretical models of the solar tachocline and 
the radiative interior currently available, a new class of helio- 
seismic instruments has to be developed, in particular to study 
the deep solar interior. Such efforts have started in Europe (e.g. 
the DYNAMICS project, Turck-Chieze et al. 2005). 
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